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ABSTRACT 

We describe and test a new method for creating initial conditions for cosmological N- 
body dark matter simulations based on second-order Lagrangian perturbation theory 
(21pt). The method can be applied to multi-mass particle distributions making it 
suitable for creating resimulation, or 'zoom' initial conditions. By testing against an 
analytic solution we show that the method works well for a spherically symmetric 
perturbation with radial features ranging over more than three orders of magnitude in 
linear scale and eleven orders of magnitude in particle mass. We apply the method and 
repeat resimulations of the rapid formation of a high mass halo at redshift ~ 50 and 
the formation of a Milky- Way mass dark matter halo at redshift zero. In both cases 
we find that many properties of the final halo show a much smaller sensitivity to the 
start redshift with the 21pt initial conditions, than simulations started from Zel'dovich 
initial conditions. For spherical overdense regions structure formation is erroneously 
delayed in simulations starting from Zel'dovich initial conditions, and we demonstrate 
for the case of the formation the high redshift halo that this delay can be accounted 
for using the spherical collapse model. In addition to being more accurate, 21pt initial 
conditions allow simulations to start later, saving computer time. 

Key words: cosmology: theory - methods: N-body simulations 



1 INTRODUCTION 

Computer simulations of cosmological structure formation 
have been crucial to understanding structure formation par- 
ticularly in the non -linear regime. Early simulations e.g. 
lAarseth et all (|l979l ) used only around a thousand particles 
to model a large region of the Universe, while recent simula- 
tions of structure formation have model led over a billion par- 
ticles in just a sin gle virialised object (|Springel et al J 1200a ; 
IStadel et aHl2009h . Since the advent of the Cold Dar k Mat- 
ter (CDM) model jPeeblej [l982l ; iDavis et al.lll985h . most 
computational effort has been expended modelling struc- 
ture formation in CDM universes. Early work in the 1980s 
focused on 'cosmological' simulations where a representa- 
tive region of the uni verse is modelled, suitable for studying 
large-scale structure (|Davis et al.l[l9 85). The starting point 
for these CDM simulations, the initial conditions, are Gaus- 
sian random fields. The numerical techniques needed to cre- 
ate the initial conditions for cosmological simulations were 
develo ped in the 1980s and are described in lEfstathiou et all 
{ 19851 ). 

As the algorithms for N-body simulations have im- 
proved, and the speed of computers increased exponen- 
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tially with time, it became feasible in the 1990s to sim- 
ulate the formation of single virialised halos in the CDM 
model with enoug h particles to be able to probe their inter- 
nal s t ructure (e.g. iNavarro et al.l 1 19961 . Il997l ; IChigna et al.l 
1 19981 ; M 

oorc et al.l Il999h .~ These more focused simulations 
required new met hods for generating the in itial conditions. 
The algorithm of iHoffman fe Ribald |l991) for setting up 
constrained Gaussian random fields was one method which 
could be applied to setting up initial conditions for a rare 
peak where a halo would be expected to form. The essence 
of this technique is that it allows selection of a region based 
on the properties of the linear density field. However, the ob- 
jects we actually observe in the Universe are the end prod- 
ucts of non-linear structure formation and it is desirable, if 
we want to understand how the structure we see formed, to 
be able to select objects on the basis of their final properties. 

This requirement led to an alternative method for 
producing initial conditions based firs t on select i ng ob - 
jects from a completed simulation e.g. iKatz et all (|l994h . 
(Navarro fc White! (|l994l ). The initial conditions for the first 
or pare nt simulation, we r e crea ted using the methods out- 
lined bv lEfstathiou et al.l (|l985h . The density field is created 
out of a superposition of plane waves with random phases. 
Having selected an object at the final (or any intermedi- 
ate) time from the parent simulation a fresh set of initial 
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conditions with higher numerical resolution in the region of 
interest, which we will call 'resimulation' initial conditions 
(also called 'zoom' initial conditions) can be made. Particles 
of different masses are laid down to approximate a uniform 
mass distribution, with the smallest mass particles being 
concentrated in the region from which the object forms. The 
new initial conditions are made by recreating the the same 
plane waves as were present in parent simulation together 
with the addition of new shorter wavelength power. 

An alternative technique for creating resimu lation ini- 
tial conditions was dev i sed bylBertschingerl J200ll ) based on 
earlier work bv lSalmonl (|l996j i andlPenl |l997lT where a Gaus- 
sian random field with a particular power spectrum is cre- 
ated starting from a white noise field. Recently parallel code 
versions using this method has been developed to generate 
very large initial conditions (|Prunet et al et al.l 

2009). A feature common to both these techniques, to date, 
is that the displacements and velocities a re set using the 
Zel'dovich approximation (IZel'dovichl 1 19701 ) . where the dis- 
placements scale linearly with the linear growth factor. 

It has long been known that simulations starting from 
Zel'dovich initial conditions exhibit transients and that care 
must be taken choosing a sufficiently high start redshift so 
that these transients ca n decay to a negligible amplitude (e.g 
lEfstathiou et al.1ll985T ). A study of the behav iour of tran- 
sients using Lagrangian perturbation theory bv lScoccimarrol 
(1998), showed that for initial conditions based on second- 
order Lagrangian perturbation theory, the transients are 
both smaller and decay more rapidly than first- order, or 
Zel'dovich, initial conditions. IScoccimarrol (|l998l ) gave a 
practical method for implementing second-order Lagrangian 
initial conditions. The method has been implement ed in 
codes for c reating cosmo l ogical initial conditions by ISirkol 
(|2005h and ICrocce et all (|2006l ). These codes are suitable 
for making initial conditions for cosmological simulations, 
targeted at large-scale structure, but they do not allow the 
creation of resimulation initial conditions, the focus of much 
current work on structure formation. 

In this paper we describe a new method for creat- 
ing second-order Lagrangian initial conditions which can be 
used to make resimulation initial conditions. The paper is 
organised as follows: in Section [2] we introduce 21pt theory 
initial conditions and motivate their advantages for study- 
ing structure in the non-linear regime by applying them to 
the spherical top-hat model; in Section [3] we describe how 
Zel'dovich resimulation initial conditions are made, and the 
new method for creating 21pt initial conditions, and test the 
method against an analytic solution for a spherically sym- 
metric perturbation; in Section U we apply the method and 
analyse the formation of a dark matter halo at high and low 
redshift for varying starting redshifts for both Zel'dovich 
and 21pt initial conditions; in Section \5\ we summarise and 
discuss the main results; in the Appendix we evaluate two 
alternate interpolation schemes which have been used in the 
process of making resimulation initial conditions. 



THE MOTIVATION FOR USING SECOND 
ORDER LAGRANGIAN PERTURBATION 
THEORY INITIAL CONDITIONS 



In this Section we give a brief summary of Lagrangian 
perturbation theory (Subsection 12. 1[) , concentrating on the 
second-order approximation, and illustrate the difference 
they make for the spherical top-hat model in Subsection l2.2l 
The method we describe could be developed further to cre- 
ate third order Lagrangian perturbation initial conditions. 
However. IScoccimarrol (1993) concluded that the gains from 
going from 2nd to 3rd order were relatively modest while 
the complexity of computing the third order terms is con- 
siderable. 



2.1 Summary of second-order Lagrangian 
perturbation theory 

We give only a bare-bones account of second-order pertur- 
bation theory here. Detailed discussions of Lagrangian per- 
tur bation theory i n general can be found in the literature 
fe.g|Bucherljll994l ; iBuchert et allll994l ; iBouchet et"ai1ll995l ; 
IScocc imarro 1998). We fo llow the notation used in Appendix 
Dl of lScoccimarrol (|l998l 1. 

In Lagrangian perturbation theory the Eulerian final 
comoving positions x are related to the initial positions q 
via a displacement field, vE'(q): 



x = q + <J>. 



(1) 



It can be shown for second-order Lagrangian perturba- 
tion theory (21pt) that the displacement field is given, in 
terms of two potentials t/r 1 ' and (f>^ 2 \ by: 

x = q-D 1 V (? <?(> (1) + D 2 V q (2) , (2) 

where D\ is the linear growth factor, and D 2 the second 
order growth factor, which is given by D 2 ~ -3D 2 /7. The 
subscripts q refer to partial derivatives with respect to the 
Lagrangian coordinates q. Similarly, the particle comoving 
velocities v are given to second order by: 



v = -D 1 f 1 HV q 



+ D 2 f 2 HV q 



kW 



(3) 



where H is the Hubble constant, the fi = dhiDi/dlna 
(a is the expansion factor). For flat models with a non- 
zero cosmological cons tant, the following relations apply 
fi w fi 5/9 , f 2 « 2fi 6/n j|Bouchet et al.ll 19951 ). where fi is the 
matter density. In practice, these approximations for fi,f 2 
and D 2 are extremely accurate when applied to making most 
actual ACDM initial conditions, as the start redshift is high 
enough that fi is very close to unity. 

The potentials <jP^ and (f>^ are obtained by solving a 
pair of Poisson equations: 

V 2 (1) (q) = 5 (1) (q), (4) 

where 8 1 \q) is the linear overdensity, and 

V 2 ^ (2) (q) = 5 (2) (q). (5) 

We will refer to the term (5"'(q), as the 'second-order 
overdensity'. In fact, this quantity is related to the linear 
overdensity field by the following quadratic expression: 



:«o=£ 



(q) 



kW, 



q)-[<#j(q) 



(6) 



where we use 4>,ij = d 2 4> / dqidqj for short. 

The most common technique used in making cosmo- 
logical i nitial conditions, to date, is the Zel'dovich approx- 
imation l|Zel' dovichlll970h . Zel'dovich initial conditions can 
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equivalently be called first order Lagrangian initial condi- 
tions, and are given by ignoring the terms with (fy 2 ^ in equa- 
tions (J2J) and ((3}. A variant of Zel' d ovich initial conditions 
was advocated by lEfstathiou et al, | i|l985l) . In this case the 
velocities are not set proportional to the displacements, but 
proportional to the gravitational force evaluated for the dis- 
placed particle distribution. This method for s etting the ve- 
locitie s gives weaker transients, as discussed in IScoccimarrol 
i|l99Sl) . although the scaling of the transients with redshift is 
unaffected. In this paper we use the term Zel'dovich initial 
conditions to mean the velocities are proportional to the dis- 
placements. We note, and we will exploit the fact later, that 
the individual terms on the right-hand side of equation ((6} 
are the derivatives of the first order Zel'dovich displacements 
which are linear in the potential <jP~' . 



2.2 The Spherical top hat model. 

The spherical top-hat model (Peebles 1980) provides a con- 
venient way to illustrate the relative merits of Zel'dovich 
and 21pt initial conditions. Apart from the significant ad- 
vantage of mathematical convenience, there are other rea- 
sons for considering this special case: 

(i) The first collapsed structures in a cosmological sim- 
ulation will form from the regions with the highest linear 
overdensity. In the limit of high peak height, the regions 
around the maxi ma in a smoothed lin ear overdensity field 
become spherical l|Bardeen et al.| [l986). As we will see, it is 
the first structures which are most sensitive to whether the 
initial conditions are Zel'dovich or 21pt. 

(ii) It is easy to show that the second-order overdensity, 
S^ 2 \ is related to the first order overdensity, through 
the inequality: 5^ < §(5^) 2 , with equality occurring 
for isotropic compression or expansion. For a fixed linear 
overdensity, the second order overdensity is maximised for 
isotropic compression, so the spherical top-hat is the limit- 
ing case where the affects of second-order initial conditions 
are maximised. 

For convenience we will consider, for now, the top-hat 
in an Einstein-de Sitter universe. In this case the spherical 
top-hat solution is most elegantly expressed in a parametric 
form as follows. For an overdense spherical top-hat which 
begins expanding at t = and reaches a maximum physical 
radius of expansion, Ro, and collapses back to zero radius 
at time, t c , the parametric solution is: 



Ro n 



cos rf) ; t 



tc_ 

2tt 



(V - sin 77) , 



(7) 



where < 77 < 2-k. These equations can be derived by solving 
the equation of motion for the top-hat: 



■k 2 RI 1 
Itl r 2 



(8) 



and applying the appropriate boundary conditions. 

If we were to solve the spherical top-hat problem in- 
stead by the techniques used in cosmological numerical sim- 
ulations, we would begin by creating initial conditions con- 
sisting of a set of particles with given masses, positions and 
velocities, at a time, t, soon after the start of the expansion. 
The solution would be obtained by evaluating the gravita- 
tional interactions between the particles and integrating the 



equations of motion of for all the particles forward in time. 
For a real simulation the accuracy of the subsequent solu- 
tion would depend not only on the accuracy of the initial 
conditions, but also on numerical aspects such as the time- 
integration scheme and the force accuracy. For our purposes, 
we will assume that these numerical sources of error can be 
overcome, and will only consider the effects of using either 
Zel'dovich or 21pt initial conditions on the subsequent evo- 
lution of the top-hat. 

The Zel'dovich and 21pt initial conditions needed to 
start a cosmological simulation of a top-hat are most easily 
obtained using the rather less elegant power-series represen- 
tation of the spherical top-hat solution. The power-series in 
t corresponding to equation ([7]) is: 



Ro 
4 



/12ttA 3 1 /127rf\3 3 /127rf 
[TTj ~ 20 ~ 2800 \~TT 



and the corresponding velocity, r is: 



(9) 



2tt.Ro 



[TTJ ~ 10 \~)~ 



2800 



127ri 



(10) 

The leading order term in equation ([9]) represents the 
expansion of a perturbation with mean cosmic density which 
grows in time in the same way as the scale factor for an 
Einstein-de Sitter universe. The second to leading term 
gives the linear growth of an overdense spherical grow- 
ing mode perturbation and it is easy to verify from equa- 
tion ([9]) that the linear overdensity at t = t c , when the 
non-linear solution collapses to a point, is given by S c — 
3(12 7r) 2/3 /20) ~ 1.686, fami liar from Press-Schechter the- 
ory (|Press fc Schechteij[l973 ). 

We obtain the Zel'dovich initial condition for the spher- 
ical top-hat by taking just the first two leading terms in r 
and r in equations @ and (|10p . The equivalent 21pt initial 
conditions correspond to taking the first three terms instead. 
We can follow the top-hat solution starting from these ini- 
tial conditions, up until the point when the sphere collapses 
to zero radius, by integrating the equation of motion equa- 
tion JS} forward in time from time, t, using the values of r 
and r at time t as the boundary conditions. 

In Figure [1] the black solid line shows the parametric 
top-hat solution given by equation 0. The red solid curve 
shows a spherical top-hat model integrated forward from 
the arbitrarily chosen time, t = t c /27, starting from the 
Zel'dovich initial condition (marked with a red dot). The 
cyan dashed curve shows the corresponding result starting 
with the 21pt initial condition at t = t c /27. The solution 
using the 21pt initial condition is very much closer to the 
analytic solution than the curve started from the Zel'dovich 
initial condition. 

The collapse time for the top-hat for both the Zel'dovich 
and 21pt initial conditions occurs later than the analytic so- 
lution. This delay, or timing error, At, can be quantified as 
follows. Multiplying equation |(HJ by r and integrating gives 
an integral of the motion: E = r 2 /2 - -K 2 Rl/2t 2 c r. For the 
analytic solution this can easily be evaluated at the radius 
of maximum expansion which gives E = —n 2 R 2 l /2t 2 . The 
corresponding values of E for the Zel'dovich and 21pt solu- 
tions can be obtained by substituting the first two or three 
terms respectively from equation ((9} and equation (|10[l into 
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Figure 1. Evolution of a top-hat spherical perturbation. The 
black curve shows the radius of a spherical top-hat perturbation as 
a function of time for an Einstein-de Sitter universe. The analytic 
solution starts expanding at t = and collapses at a time t = t c 
and has a maximum physical radius of expansion of i?o- The 
red and dashed cyan curves show the spherical top-hat solution 
obtained by integrating the equation of motion for the top-hat 
starting from t = t c /27 and respectively initialising the radius 
and its velocity with the Zel'dovich and 21pt initial condition. 
The result for the 21pt initial condition is more accurate than 
for the Zel'dovich initial condition where the time of collapse is 
significantly delayed. 



the expression for E. The expansion time, t c , of the top-hat 
over the complete cycle, by analogy with Keplerian orbits, 
is related to E by: dlni c /dln \E\ = — | for bound motion in 
an Einstein-de Sitter universe. We can generalise this result 
to fiat universes with a cosmological constant by numeri- 
cally integrating equation ((8| with an added term for the 
cosmological constant. In this case we find, by fitting to the 
results of numerical integration, an equivalent approximate 
relation: d!nt c /dln|S| = — |r2°' 23 , valid to an fractional 
accuracy of less than 1.2 percent for 0.25 < fi c < 1, where 
Q c is the matter density in units of the critical density at 
the time of collapse, provided t s <C t c so the cosmological 
constant is negligible at t s . 

Knowing dlnt c /dln \ E\ and assuming again the condi- 
tion that the starting time of the initial condition, t s , obeys, 
t s <C t c , allows the fractional timing error, At/t c , to be ob- 
tained, to leading order in t s /t c , for both types of initial 
condition, from the ratio of E for the analytic solution with 
respectively with the value of E for Zel'dovich or 21pt initial 
conditions. 

Applying these results we find after some algebra that 
the timing error for Zel'dovich initial conditions to lowest 
order is: 



At 

77 



35e 

10 



2/3 

- 1 n: 



Similarly, the 21pt initial condition fractional timing error is 
given by: 



At 23S 2 C /0 4/3 



tc 



210 \t 



- n: 



(12) 



Alternatively, these results can be expressed in terms 
of the start redshift of a numerical simulation, z s . Relative 
to an analytic spherical top-hat with collapse occurring at 
redshift, z c , where z c <C z 3 , the fractional timing error is 
given for Zel'dovich initial conditions by: 



At 

77 



0.51 fl + z 



fio- 36 \l + z, 
and for 21pt initial conditions it is: 

0.31 fl + z, 



At 
77 



ftg- 49 \l + z 



(13) 



(14) 



(11) 



where we have made use another approximate relation: 
t 2/3 (l + z) oc f2 013 relating redshift, z, to the age of the 
universe, t and the value of Q all at time, t, which is accu- 
rate to better than one percent in the range 0.25 < fi < 1 
for a flat universe with a cosmological constant. 

Ignoring the weak fi dependence, the timing error for 
the 21pt initial condition shows a quadratic scaling with the 
ratio of the expansion factor at the start to that at the time 
of collapse, while for Zel'dovich initial conditions, the scal- 
ing is just linear. These results for spherical collapse are 
consistent with the mor e general behaviou r of the decay 
of transients discussed in IScoccimarrd (1998). The absolute 
value of the timing error, At, for a fixed start redshift, ac- 
tually decreases with increasing collapse time for 21pt initial 
conditions. For Zel'dovich initial conditions, the timing er- 
ror grows slowly with increasing collapse time. However, by 
starting at sufficiently high redshift, the timing error for the 
Zel'dovich initial conditions can be made small. It is common 
practice to start simulations from Zel'dovich initial condi- 
tions at a high redshift for this reason so that the transients 
have sufficient time to decay to a negligible amplitude. 

Equating the fractional timing errors for equations 
equation fl 13 p and equation (|14p gives an objective way of 
defining an equivalent starting redshift for Zel'dovich or 21pt 
initial conditions for a given choice of collapse redshift. To be 
conservative, it would be natural to choose the collapse red- 
shift to be the first output redshift of the simulation which is 
of scientific interest. To give a concrete example, ignoring the 
Q dependence, if the 21pt initial conditions are created a fac- 
tor of 10 in expansion before the first output, the equivalent 
Zel'dovich initial conditions would need to start at a factor 
of about 160 in expansion before the first output. Using 21pt 
initial conditions instead of Zel'dovich initial conditions can 
therefore lead to a significant reduction in computer time 
because the simulations can be started later. A later start 
has further advantages as force errors, errors due to the dis- 
creteness in the mass distribution and time integration er- 
rors all accumulate as the ratio of the final expansion factor 
to the initial expansion factor increases. All of these reasons 
make 21pt initial conditions attractive for the simulator. In 
the next section we describe the new method to make 21pt 
initial conditions. 
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3 MAKING SECOND ORDER LAGRANGIAN 
RESIMULATION INITIAL CONDITIONS 

In this Section we explain what resimulation initial condi- 
tions are and describe how they can be made. We will de- 
scribe the method currently in use to make initial conditions 
for the Virgo Consortium for Supercomputer Simulation^. 
The Virgo resimulation code to date generates Zel'dovich 
initial conditions. In Subsection 13. 1 1 we describe the process 
for making these Zel'dovich resimulation initial conditions. 
In Subsection 13.21 we describe the new method to generate 
21pt initial conditions. In Subsection l3.3l we demonstrate by 
comparing with an analytic solution that the method works 
in practice. 



3.1 Setting up Zel'dovich resimulation initial 
conditions 

The current method u sed by the Virgo code is based on the 
ideas first described in lNavarro fc White! (1 19941 ). The actual 
implementation has changed considerably over time with in- 
creasi n gly up-to-date accou nts b eing given in | Power et all 
l|2003t ). lNavarro et alj (|2004l ) and lSpringef et al. I < |2008h . 

The need for resimulation initial conditions stems from 
the desire to repeat a particular simulation of structure for- 
mation with (usually) higher mass resolution in one or more 
regions than was present in the original, or parent, simu- 
lation. The parent simulation may have been run starting 
from cosmological initial conditions or it may be a resim- 
ulation itself. The following list gives the steps needed to 
create a set of resimulation initial conditions from a parent 
simulation. 

(i) Identify the region(s) of interest in the parent simula- 
tion where higher mass resolution is required. 

(ii) Create a force-free multi-mass particle distribution 
which adequately samples the high resolution region(s) of 
interest identified in (i) with low mass particles. This par- 
ticle distribution needs to be a good approximation to an 
isotropic and homogeneous mass distribution. The initial 
density fluctuations are created by perturbing the positions 
and assigning velocities to each of the particles. The region 
exterior to the high resolution region can often be sampled 
more crudely than the parent simulation, using larger mass 
particles, as its only purpose is to ensure that the gravita- 
tional tidal field in the high resolution region is accurately 
reproduced. Typically the masses in the exterior region scale 
roughly as the cube of the distance from the high resolution 
region. We will refer to this special particle distribution from 
now on as the 'particle load'. 

(iii) Recreate the displacement and velocity fields for the 
parent simulation and apply the displacements and assign 
the velocities to the particle load created in (ii). 

(iv) Add extra small-scale fluctuations to the region of 
interest only, to account for the fact that the mass distribu- 
tion in the resimulation is now more finely sampled than the 
parent simulation and therefore has a higher spatial Nyquist 
frequency. For Zel'dovich initial conditions this high spatial 
frequency power can be added linearly to the low frequency 
power copied from the parent simulation. 
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In this paper we will not be concerned with steps (i) 
and (ii), which are common to resimulations in general, but 
will concentrate on adapting steps (iii) and (iv) to allow the 
creation of 21pt initial conditions. 

The displacement and velocity fields needed to perturb 
the particle load can in practice be efficiently created using 
Fourier methods. Adopting periodic boundary conditions 
when setting up the perturbations forces the Fourier rep- 
resentation of the density, displacement and velocity fields 
to be composed of discrete modes in Fourier space. The am- 
plitudes and phases of these modes are set as described in 
lEfstathiou et al.l (119851 ) to create a Gaussian random field 
with the appropriate power spectrum in Fourier space. Start- 
ing from the Fourier representation of the density field it is 
straightforward to use discrete Fourier transforms to gen- 
erate each of the components of the displacement field in 
real space. The discrete Fourier transform gives the displace- 
ments on a cubic lattice. 

In practice, because computer memory is limited, more 
than one Fourier grid is required to be able to complete 
step (iv). These grids have different physical sizes, with the 
largest grid covering the entire periodic cubic volume, as is 
the case when making cosmological initial conditions. The 
additional grid(s) are given a smaller physical size, which 
means even with limited computer memory, that they can 
generate a displacement field with power at higher wavenum- 
bers than the largest grid, but only over a restricted part of 
the simulation volume. By nesting sufficient Fourier grids 
about the high resolution region it is possible to generate 
power all they way down to the particle Nyquist frequency 
in the high resolution region. A description of the creation of 
the initial conditions for the Aquarius halos, which required 
two Fourier grids, is given in lSpringel et all (|2008T l. 

For Zel'dovich initial conditions, the displacement field 
at any location is formed by coadding the displacement con- 
tributions from the one or more of these overlapping cubic 
Fourier grids, and the corresponding velocity field is simply 
proportional to the displacement field. For each Fourier grid 
only a subset of the possible modes are assigned non-zero 
values. In k-space these 'populated' modes are split between 
the different grids in a nested fashion around the origin. Each 
grid populates a disjoint region of k-space concentric on the 
origin, and taken together the modes from all the grids fill a 
spherical region in k-space around the origin. The smallest 
physical grid, which typically fits around just the high res- 
olution region, is used to represent the highest wavenumber 
modes in k-space. The highest wavenumber possible is lim- 
ited by the particle Nyquist frequency which is determined 
by the interparticle separation. 

A displacement field is generated on each Fourier grid. 
The complete displacement field is calculated by adding the 
the contributions from the various grids together. To do this 
it is necessary to define a set of common positions at which 
to co-add the displacements. The most natural common po- 
sitions to use are those of the particle load itself, as the 
displacements need to be known at these locations to gen- 
erate the Zel'dovich initial conditions. The displacements at 
the particle positions are computed by interpolating the val- 
ues from the nearest Fourier grid points. A discussion of the 
interpolation methods currently in use and their associated 
errors is given for completeness in the Appendix. The inter- 
polation errors only become significant for fluctuations with 
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a wavelength comparable to the Nyquist wavelength of the 
Fourier grid. 

For the high resolution region a second scale, the par- 
ticle Nyquist wavelength also becomes important. Even in 
the absence of interpolation errors, the discreteness of the 
mass distribution affects the growth of linear fluctuations 
with wavelengths close to the particle Nyquist frequency 
(I Joyce et al.ll2005T ). We will not be concerned about either 
of these sources of error in the main part of the paper. 



3.2 Generating 21pt resimulation initial conditions 

In Subsection 13.11 we outlined a method for making 
Zel'dovich initial conditions for resimulations. In this Sub- 
section we will take this method for granted, and concen- 
trate on how to augment the Zel'dovich initial conditions to 
produce 21pt initial conditions. This is achieved, as can be 
seen from equations ((2| and ((3]), by computing the com- 



ponents of the field V q 



t(2) 



The second-order corrections to 



the Zel'dovich initial conditions for both the positions and 
the velocities are proportional to this quantity. 

The second-order potential, V q (j>^ 2 \ is obtained by solv- 
ing the Poisson equation ([5]), where the source term is the 
second-order density field, 5^ 2 \ given by equation {B). The 
right-hand side of this equation is a non-linear combination 
of six distinct terms, each of which are second derivatives of 
the first order potential, <jP^ . This same first order potential 
is used to generate the Zel'dovich displ acements and veloci- 
ties. For this reason, as pointed out by Sco ccimarrol fl998), 
the calculation of the second-order overdensity field, fits in 
very well with the generation of Zel'dovich resimulation ini- 
tial conditions, from a computational point of view. The six 
second derivatives can be calculated using Fourier methods 
in the same way as is used to compute the first derivatives 
of <j>^ needed for the Zel'dovich approximation. 

For initial conditions that can be made using a single 
Fourier grid, the six derivatives can be combined to give 
the value of at each Fourier grid point. The Poisson 
equation (JSJ, can then be solved immediately usi ng stan- 
dard Fourier methods <|Hocknev fc Eastwood! Il988h to give 
the components of V ' q (j>^ at the grid points. Choosing a 
suitable cubic grid arrangement for the particle load en- 
sures that the Lagrangian positions of the particles coincide 
exactly with Fourier grid points. This means the particle 
21pt displacements and velocities can be evaluated without 
the need for any spatial interpolation. This is the app roach 
adopted by the 21p t codes made available ISirkd |2005l ) and 
ICrocce et all l|2006l ). 

The approach just described does not work when more 
than one Fourier grid is used to generate the Zel'dovich 
initial conditions, because <5' 2 ' is a non-linear function of 
the first-order potential, <jP^ , unlike the Zel'dovich displace- 
ments and velocities which are linear functions. The solution 
to this problem is simply to compute all six of the deriva- 
tives, which are linear functions of <fr \ for each Fourier grid, 
and co-add them. Then once the summation is complete for 
all six derivatives for a given particle, compute 8^ using 
equation ([6|. 

As with the Zel'dovich displacements, the values of the 
six derivatives need to be interpolated to a common set of 
positions before they can be co-added. It proves convenient 
again to use the particle load to define these common po- 



sitions. Computationally this entails storing six additional 
quantities per particle in computer memory. We use the QI 
interpolation scheme, described in the Appendix, to inter- 
polate the derivatives to the positions given by the particle 
load, as this scheme helps to minimise the interpolation er- 
rors, given limited computer memory. We have implemented 
the algorithm to compute <5' 2 ' for each particle in the par- 
ticle load in a Fortran/C code called IC_2LPT_GEN. This 
code also computes the Zel'dovich displacements and veloc- 
ities for each particle at the same time. 

The final stage to make the 21pt resimulation initial con- 
ditions is to solve the Poisson equation (JSJ). At this point it 
is no longer trivial to solve the equation by Fourier methods 
because the second-order overdensity field is known only at 
the positions of the particle load, which do not in general 
form a regular grid. However it is relatively simple to solve 
the Poisson equation using an N-body code Poisson solver. 
It is natural, given that the purpose of creating the initial 
conditions in the first place is to integrate them forward in 
time with an N-body code, to use the same N-body code 
to solve both equation (0 and to integrate the equations of 
motion. 

As mentioned in Subsection 13.11 in point (ii), the La- 
grangian positions of the particle making up the particle 
load are special. The particles are distributed so that the 
mass density field is an excellent approximation to a uniform 
mass density, on scales larger than the interparticle separa- 
tion. Furthermore the gravitational force on any particle due 
to the others cancels out to a very good approximation. We 
can exploit the properties of the particle load in a simple 
way to reconstruct the second-order density field. To do this 
we create a '21pt particle load' by taking the positions from 
the particle load, but replacing the particle masses with a 
'21pt' mass, given by: 



T712LPT 



= m(l + A<5 (2) ), 



(15) 



where A is a constant numerical factor, whose value for now 
we will take to be unity. 

With minimal modification, an N-body code can then 
be used to compute the gravitational force due to the 21pt 
particle load acting on itself. By construction this amounts 
to solving equation (5), with the resulting gravitational 



acceleration giving, V, 



at the positions given by the 



particle load. Provided the N-body code also reads in the 
Zel'dovich velocities for the particle load, the 21pt initial con- 
dition positions and velocities can be evaluated using equa- 
tions ((2} and <(3j respectively. Discarding the 21pt masses, 
and replacing them with the true particle masses completes 
the creation of the 21pt resimulation initial conditions. The 
N-body code can then be used in the normal way to integrate 
the equations of motion for all of the particles. 

The amplitude of the second-order overdensity field, 
and similarly S7 q (j>^ depends on the redshift of the initial 
conditions. By introducing the scaling factor A into the '21pt' 
masses, and rescaling the resulting V,^' 2 ' by 1/A to com- 
pensate, it is possible to make the calculation performed by 
the N-body code independent of the requested redshift of the 
initial conditions. It is desirable to implement this scaling if 
one wants to calculate V q 4>^ accurately for high redshift 
initial conditions, as otherwise the inevitable force errors 
from the N-body force calculation will begin to dominate 
over the small contribution from the second-order overden- 
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sity field. The precise value of A is unimportant, provided it 
is large enough. A suitable choice is given by setting A to be 
as large as possible while avoiding any of the '21pt' masses 
becoming non-positive. 

We have implemented the final step to make the 21pt 
initial conditions in practice with the P-Gadget3 code, writ- 
ten by Volker Springel. P-Gadget3 is an N-body/SPH code 
writte n in C. This c ode is similar to the public Gadget-2 
code jSpringej [20051 ), and has been used, for example, for 
the Aquarius project (|Springel et alj 12008). Operationally 
the Fortran/C code, IC_2LPT_GEN, on completion outputs 
a special format initial condition file(s) for the P-Gadget3 
code to read in. The file is special in that the particle po- 
sitions are the Lagrangian positions, i.e. the particle load 
positions, the velocities are the Zel'dovich velocities, and an 
additional block of data is provided which gives every par- 
ticle a second mass - the '21pt-mass' described above. This 
required a modification to P-Gadget3 so that it first can 
recognise the special file format required by the 21pt initial 
conditions through a new flag stored the file header. Hav- 
ing read in the initial conditions, P-Gadget3 first does a 
force calculation using the 21pt masses for each particle as 
described above, and once the 21pt initial conditions have 
been created, integrates the equations of motion forward in 
time as usual. 

The linear density field created by the r esimulation 
metho d outlined above is bandwidth limited (see lPress et al.l 
(f992) for a discussion). The second-order overdensity field 
is similarly bandwidth limited but because equation ([6]) is 
quadratic, the second-order overdensity has twice as large a 
bandwidth in each physical dimension. If the second-order 
overdensity field is known at the Lagrangian particle posi- 
tions then it is only fully sampled if the linear density is 
limited to a wavelength which is twice that of the particle 
Nyquist wavelength (this approach is implemented in the 
ICrocce et al.l (|2006l ) code). It is common practice in mak- 
ing initial conditions to use waves for the linear density 
field down to the Nyquist wavelength itself. In this case the 
second-order overdensity field at the positions of the parti- 
cles undersamples the second-order density field. However, 
as mentioned earlier, the growth rate of fluctuations of very 
short wavelengths is also compromised by the discreteness of 
the mass distribution in any case whether using Zel'dovich 
or 21pt initial conditions so we will ignore this issue, but it 
remains the case that it is always advisable to do resolution 
tests to establish the limitations of the simulations. 



3.3 Testing the 21pt resimulation initial 
conditions 

We can test limited aspects o f the IC_2LPT_GE N code 
against the code made public bv lCrocce etafl (|2006f ) for the 
cases of both Zel'dovich and 21pt cosmological initial condi- 
tions, where the second-order overdensity field can be com- 
puted using a single Fourier grid. The two codes were writ- 
ten i ndependently alth ough they do share som e subroutines 
from | Press et al.l l|l992i ) . We find that when the lCrocce et al.l 
(2006) code is modified so as to set up precisely the same 
density field (i.e. exactly the same modes, with the same 
phases and amplitudes for each mode, and the same shape 
linear input power spectrum, the same power spectrum nor- 



malisation and the same cosmological parameters) then the 
output of the two codes agrees extremely well. 

Testing the 21pt resimulation initial conditions gener- 
ated using the method described in the previous section for 
a general case is problematic as it is non-trivial to compute 
the second order terms by alternative means. However, it is 
relatively simple to test the code for a spherically symmetric 
perturbation for which one can calculate the properties of 
the initial conditions analytically. The IC_2LPT_GEN code 
needs only minimal alteration to be able to create a spher- 
ical perturbation. The only change needed is to replace the 
part of the code which generates the amplitudes and random 
phases for each Fourier mode, as required for generating ran- 
dom Gaussian fields, with a module that instead computes 
the appropriate amplitude and phase for a specified spherical 
perturbation with a centre of symmetry at a given location. 
In other respects the code is unaltered, which means the 
spherical test remains a powerful test of the overall method 
and the particular numerical implementation. 

The test initial condition most closely resembles what 
would be expected around a high overdensity peak and is 
similar to the initial conditions in Section [4] for t he highest 
resolu tion initial cond i tions used in the papers iGao et all 
(2005) and lReed et~aH j|2005l ). These two papers were based 
on a simulation of a relatively massive and very rare high 
redshift halo, a suitable location for finding one of the first 
metal free stars. 

For our test we add a spherically symmetric pertur- 
bation, to a homogeneous and isotropic Einstein de-Sitter 
universe. We adopt a simple analytic form for which the 
second-order overdensity and the 21pt displacements and ve- 
locities can be computed analytically. The perturbation was 
formed by summing n component terms to give the linear 
overdensity field as follows: 

n 

S W (r) = XX W (3-r 2 /a, 2 )exp(-r 2 /2a|), (16) 
j=i 

where r is the radius, 5^ is the amplitude of the linear 
overdensity, and Oi the linear size, for the ith component of 
the perturbation. The prefactor to the exponential is chosen 
so that: r 2 5^'dr = 0, meaning that the net mass of the 
perturbation is zero. 

The choice of a radial Gaussian cut-off ensures that 
the spherical perturbation is well localised which is a nec- 
essary requirement: the use of periodic boundary conditions 
in making the resimulation initial conditions violates spher- 
ical symmetry, and means the choice of the scale lengths <jj 
has to be restricted or the resulting perturbation will not 
be close to being spherically symmetric. We have chosen the 
scale lengths, <Ji, so that each component of the perturba- 
tion is essentially generated by a corresponding Fourier grid, 
with the value of a% being no greater than about 10 percent 
of the side length of the grid and each grid is centred on 
the centre of spherical symmetry. These requirements en- 
sure that each component of the perturbation is close to, 
although not exactly, spherical. 

All that is needed to build the perturbation with the ini- 
tial conditions code is the equivalent Fourier representation 
of the real space perturbation equation (|16[) . The Fourier 
representation, <5(k) for a perturbation with overdensity 5(x) 
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is given by: 

5(k) = j S(x) exp(ik ■ x)d 3 x, (17) 

where k is the wavevector. Equation (|16[) has the following 
Fourier representation as a function of the magnitude of the 
wavevector k: 

n 

S(k) = (2tt) 3/2 E SfVff? exp(-fc 2 a 2 /2). (18) 
The solution of equation (j4| is: 

n 

<P W (r) = £ «J V exp(-r 2 /2a 2 ). (19) 

i = l 

For a general spherical perturbation, the source term for the 
Poisson equation (0 simplifies in spherical polar coordinates 
to: 

which can be integrated to give the second order displace- 
ment: 

€ = i^) S , (21) 
dr r \ dr J ' 

where the constant of integration has been set to zero to 
avoid a source term at the coordinate origin. 

For the spherical perturbation given by equation (|16[) . 
the potential is given by: 

n n 

<^ (2) M - EE^^4« P (V/4), (22) 
i=i j=i 

where 2/cr 2 ,- = 1/cr 2 + 1/cr 2 . 

The second-order overdensity, given by equation © is: 

S (2 \r) = £E W (3 ^ Sf) exp(-r 2 /4). (23) 

For the test, a set of initial conditions was made sum- 
ming five terms in Eq. (|16|l . The simulation volume was 
chosen to be periodic within a cubic volume of side length 
100 units and the centre of symmetry was placed in the 
middle of the simulation volume. Five Fourier grids were 
used also centred on the centre of the volume. The param- 
eters chosen are given in Table [T] The Table gives the val- 
ues of Sj and <n for the five components making up the 
spherical perturbation. For the purposes of the test the ac- 
tual value of the overdensity is unimportant. We arbitrar- 
ily chose the central overdensity to be unity. This is an or- 
der of magnitude higher than what would typically occur 
in a realistic simulation. The quantity L gr id gives the phys- 
ical dimension of each Fourier grid. The first grid, which 
covers the entire cubic periodic domain, has a dimension 
of 100 units on a side. The parameter JV gr id is the size 
of the Fourier transform used in the computations. The 
parameters K m i B and if max define which Fourier modes 
are used to make the spherical perturbation. All modes 
are set to zero except for the modes in each grid with 
wavenumber components, 2n/L gT i,i(k x ,ky,k z ), which obey 
-Kmin < max(|fc x |, \k v \, | fe z |) < -Kmax- These limits apply to 
grids 1, 2, 3, 4 and grid 5, except for the high-k limit, which 
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Figure 2. Comparison between quantities generated by the rcs- 
imulation codes for a spherical perturbation and the exact ana- 
lytical result. The top panel shows the Zel'dovich velocities (uW) 
are accurately reproduced by the IC_2LPT_GEN code. The mid- 
dle panel shows the second-order overdensity (<5' 2 )), output by 
the IC_2LPT_GEN code. This quantity is needed as an input to 
the P-Gadget3 code to generate the second order Lagrangian dis- 
placements and velocities. The bottom panel shows the second 
order Lagrangian velocities computed by the P-Gadget3 code us- 
ing the values of the second-order overdensity computed for each 
individual particle. The arrows mark the side lengths of the five 
periodic Fourier grids that are used to generate the perturbation 
in the rcsimulation code. The dotted vertical lines mark the short 
wavelength cut-off for each of the Fourier grids. 



is replaced by the condition v/fcl + ky + fc| < T^max to en- 
sure a spherical cut-off in k-space to minimise small-scale 
anisotropy. 

Following stage (ii) of the process described in Subsec- 
tion [371] a special particle load was created for the test con- 
sisting of concentric cubic shells of particles centred on the 
centre of the volume. In detail, the cubic volume of side 
length 100 was divided into n 3 cubic cells. A particle was 
placed at the centre of just those cells on the surface of the 
cube (a total of n 3 — (n — 2) 3 cells). The mass of the particle 
was given by the cell volume times the critical density. The 
remaining (n — 2) 3 cells define a smaller nested cube of side 
length 100 x (n — 2)/n. This smaller cubic volume was again 
divided into n 3 cells and the process of particle creation re- 
peated. The process was repeated 360 times, taking n = 51. 
On the final iteration all n 3 cells were populated, filling the 
hole in the middle. Provided the spherical distribution is 
centred on the centre of the cube this arrangement of parti- 
cles allows a spherical perturbation covering a wide range of 
scales to well sampled everywhere. The particle load has cu- 
bic rather than spherical symmetry, so this arrangement will 
generate some asphericity, but for sufficiently large n, this 
becomes small. Similar results were obtained with a smaller 
value, n — 23. 

Figure [2] shows the results of the test. The top panel 
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0.1 


15 


100 


512 


1 


80 


2 


0.15 


2 


15 


512 


12 


80 


3 


0.2 


0.2 


1.5 


512 


8 


80 


1 


0.25 


0.024 


0.18 


512 


10 


80 


5 


0.3 


0.003 


0.0225 


512 


10 


80 



Table 1. Parameters for a spherical perturbation used to test 
the method for generating 21pt initial conditions. The index i 
and quantities 6, and Oi characterise the perturbation used for 
the test as given by equation dl6l . The remaining columns give 
parameters of the Fourier grids used to realise the perturbation 
and are explained in the main text. 



shows the Zel'dovich velocities (i/ 1 ') as a function of ra- 
dius, scaled by the Hubble velocity relative to the centre. 
The black line shows the analytic solution, while the red 
dots show binned values measured from the output of the 
IC_2LPT_GEN code. The vertical arrows mark the edge 
length of the five cubic Fourier grids, while the five ver- 
tical dotted lines show the high spatial frequency cut-off 
(given by K max ) of the imposed perturbations for the cor- 
responding cubic grid. The first panel establishes that the 
IC_2LPT_GEN code is able to generate the desired linear 
perturbation. 

The second panel of Figure [2] shows in black the second- 
order overdensity of the analytic solution, while the red dots 
show the results from the output of the IC_2LPT_GEN. This 
shows that the second-order overdensity is successfully re- 
produced by the code. 

The third panel of Figure [2] shows the second-order ve- 
locity (i> )• The black curve is the analytic solution and the 
red dots show the binned second-order contribution gener- 
ated by the P-Gadget3 code, using the values of the second- 
order overdensity computed for each particle which are gen- 
erated by the IC_2LPT_GEN code. The method works well, 
the fractional error in the value of is below 1 percent 
everywhere except close to the scale of the first grid where 
spherical symmetry is necessarily broken and the perturba- 
tion has an extremely small amplitude anyway. 

In conclusion the method is able to work well for a per- 
turbation that spans a very large range of spatial scales and 
particle masses. The 5th component has a length-scale of 
just l/5000th of the 1st component and the particle masses 
associated with representing these components vary by the 
cube of this number - more than eleven orders of magnitude. 
Most practical res imulations are le ss demanding, although 
the resimulation of lGao et al.ll|2005h . which we reproduce in 
the next Section, is similar. 



4 REAL APPLICATIONS 

In this Section we remake some resimulation initial condi- 
tions with Zel'dovich and 21pt initial conditions and run 
them with the P-Gadget3 code. We investigate how the end 
point of the simulation depends on the start redshift and 
the type of initial condition. 




Z.6xl0»M o /h ■ 

6.4x1 5 M o /h - 

1.6xl0 6 M e /h 

4xlO«M B /h 
lO^/h 
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0.01 0.015 0.02 0.025 

Expansion factor 

Figure 3. The physical radius of a sphere containing a fixed 
amount of mass in a simulation of the formation of a massive halo 
at high redshift. The sphere is centred on a density maximum 
located by applying the shrinking sphere algorithm detailed in 
the text. The four curves show simulations starting from either 
Zel'dovich or 21pt initial conditions, at a start redshift of 399 or 
599. The agreement between the two 21pt initial conditions is very 
good. As expected from the spherical collapse model, discussed 
in Subsection [22] the collapse of the Zel'dovich initial conditions 
is delayed with the largest delay occurring for z s = 399. 



4.1 Structure formation at high redshift 

The initial conditions for the spherically symmetric test in 
Subsection l3.3l was intended to approximate the initial con- 
dition s in a series of resimulations presented in iGao et all 
i|2005h to investigate the formation of the first structure in 
the Universe. In that paper, a cluster was selected from a 
simulation of a cosmological volume. A resimulation of this 
cluster was followed to z = 5, and the largest halo in the 
high resolution region containing the proto-cluster was it- 
self resimulated and followed to z = 12. The largest halo in 
the high resolution region at this redshift was then resimu- 
lated and followed to z = 29. This was repeated one more 
time ending with a final resimuation (called 'R5') yielding 
a halo of mass M200 = 1.2 x 10 5 h~ l Mq (that is the mass 
within a spherical region with a mean density of 200 times 
the critical density) at redshift 49, simulated with particles 
with a mass of just O.545/i _1 M0. The original motivation 
for studying such a halo is that halos like this one are a po- 
tential sit e for the formatio n of the first generation of metal 
free stars l|Reed et al.ll2005T ). 

We have recreated the R5 initial conditions with 
both Zel'd o vich a nd 21pt resimulation initial conditions. In 
IGao et al l (|2005T ). the starting redshift was 599, and we 
have repeated a simulation starting from this redshift and 
also z s — 399, for both Zel'dovich and 21pt initial condi- 
tions. We ran each simulation with P-Gadget3 to z = 43| 
and produced outputs at expansion factors of 0.01, 0.015, 
0.017, 0.019, 0.021, 0.0225. These initial conditions remain 
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Figure 4. Same as Figure [3] except that the expansion factor 
plotted for each of the four curves has been corrected to account 
for the time delay expected in the spherical collapse model and 
given by equations J 13H and J 14t . The density inside each radius 
exceeds ten times the critical density to the right of the black 
dots. The correction is expected to work well for a spherical col- 
lapse above an overdensity of around ten (see main text). Even 
through the collapse is not perfectly spherical applying the cor- 
rection does greatly improve the match between the four curves 
particularly to the right of the black dots. The spacing of the 
outputs is insufficient to resolve features around a = 0.018 where 
shell crossing becomes important in the lower two curves. 



the most complex that have been created by the Virgo 
resimulation code, requiring seven Fourier grids. The very 
large dynamic range also required that the simulation of 
R5 itself was carried out with isolated boundary conditions 
following a spherical region of radius 1.25/i _1 Mpc which 
is much smaller than the parent simulation of side length 
479fa _1 Mpc. We adopted the same isolated region for the 
new simulations. 

The use of isolated boundary conditions gives the iso- 
lated region as a whole a velocity which depends both on 
the starting redshift and initial condition type. In order to 
compare the different versions of R5 at particular redshifts 
it is first necessary to establish a reference point or centre 
from which to measure. Using a position determined from 
th e main hal o is on e potential approach, but it was found 
in iGao et all (|2005T l that using a group finder like friends- 
of-friends on this rather extreme region results in groups 
that are highly irregular in shape, and are therefore unsuit- 
able. Instead we have found a centre for each simulation 
output usin g the shrinking sphere algorithm described in 
I Power et al.l l|200&t ) starting with a sphere that encloses the 
entire isolated region. Having found a centre, we have grown 
spheres from this point to find the radii enclosing masses of 
( 10 4 , 4 x 10 4 , 16 x 10 4 , 64 x 10 4 , 256 x 10 4 ) h~ 1 M . We checked, 
for one output redshift, by making dot plots that this centre 
can be identified visually as being in the same recognisable 
lump. 



Figure [3] shows these five radii as a function of the ex- 
pansion for the four resimulations. The radii for the two res- 
imulations starting from 21pt initial conditions at redshifts 
399 and 599 are extremely close. The simulations starting 
from Zel'dovich initial conditions show a significantly de- 
layed collapse, with the lower start redshift showing the 
largest delay. This is in accordance with what would be ex- 
pected from the predictions of the spherical collapse model 
as shown in Subsection 12.21 

The collapse occurring in the R5 simulation is only ap- 
proximate l y sph erical as evidenced by Figures 3 and 4 in 
IGao et ahl (|2005l ). It is nonetheless interesting to see whether 
the results from the spherical top-hat model do agree quan- 
titatively as well as qualitatively. To do this we first esti- 
mated an expansion factor of collapse for each of the five 
21pt z 3 — 599 curves, by fitting a spherical top-hat to each 
curve, and taking a value for the expansion factor of collapse 
from the corresponding top-hat model. We chose these par- 
ticular curves because the expected timing error for these 
21pt initial conditions should be very small, compared to 
the simulations started from Zel'dovich initial conditions, 
and so should not bias the estimate of the collapse expansion 
factor significantly. In Figure [4] the same radii are plotted 
as in Figure However the values of the expansion factors 
along the x-axis are 'corrected' using the equations (|13[) and 
(|14|) to cancel out the expected timing errors for both the 
21pt and Zel'dovich initial conditions. The timing error is 
assumed constant along each curve. Strictly speaking this 
correction applies at the moment the spherical top-hat so- 
lution collapses to a point. However inspection of Figure [T] 
together with the fact that top-hat solution is symmetric in 
time, for any initial condition, about the point of maximum 
expansion (for an Einstein-de Sitter universe, an excellent 
approximation at redshift ~ 50.), shows that the timing er- 
ror mostly arises around the time of maximum expansion 
and therefore the correction at the time of collapse, should 
be applicable once the top-hat has collapsed by more than 
~ 0.85 of the radius of maximum expansion. For the spher- 
ical top-hat this corresponds to the density exceeding the 
critical density by a factor of about ten, and this overden- 
sity is marked on each curve in Figure U by a black dot. 
To the left of the dots the density is less than ten so the 
correction is expected to be too strong, which it is. To the 
right of the dots the shift works rather well with the lines 
showing much better agreement than in Figure [3] The spac- 
ing of the outputs is too crude to follow the collapse of the 
innermost two radii very accurately, so it is not surprising 
the agreement around a — 0.018 is less good than elsewhere. 

In conclusion, the results using 21pt initial conditions 
are much less sensitive to the start redshift and it is possi- 
ble to start the simulation later. The formation of the halo 
is delayed using Zel'dovich initial conditions and the mag- 
nitude the delay agrees with what one would expect for a 
spherical top-hat collapse. 

4.2 Aquarius halos from 21pt initial conditions 

As a final example, we look at the effect of using 21pt resim- 
ulation initial conditions on the r edshift zero internal struc- 
ture of one of the Aquarius halos |Springel et al.| [2008V The 
start redshift of the Aquarius halos was 127, so even for 
a spherical top-hat collapse the timing error given by equa- 
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Figure 5. The circular velocity as a function of radius for the 
Aquarius Aq-A-5 halo at z=0 for simulations starting from 21pt 
initial conditions (shifted vertically by 30 kms -1 for clarity) and 
Zel'dovich initial conditions, starting from the redshifts indicated 
in the key. The curves for 21pt initial conditions arc much more 
consistent amongst themselves. The Aquarius Aq-A-5 halo in 
Springel et al. (2008) was started at z s = 127 from Zel'dovich 
initial conditions. In the absence of the 30 kms^ 1 vertical shift, 
the z B = 127 curve starting from Zel'dovich initial conditions lies 
between or very close to the 21pt circular velocity curves. 



tion would only be around 90 Myrs for Zel'dovich initial 
conditions, which is short compared to the age of the Uni- 
verse, so any effects on the main halo are expected to be 
small. 

We chose th e Aq-A-5 halo, with a bout 600000 particles 
within 7-200, from lSpringel et al.l l|2008l ) as a convenient halo 
to study the effects of varying the start redshift and initial 
conditions. For this study, we ran eight (re)simulations of 
Aq-A-5 to redshift zero, starting at redshifts 127, 63, 31 and 
15, for both Zel'dovich and 21pt initial conditions. We used 
the P-Gadget3 (which is needed to create the 21pt initial 
conditions) and took i dentical values f or the various numer- 
ical parameters as in ISpringel et al. (2008), and similarly 
ran the SUBFIND algorithm to identify subhalos. At this 
resolution, there are typically about 230 substructures with 
more than 31 particles in the main halo. 
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Figure 6. Cumulative subhalo abundance as a function of the 
maximum subhalo circular velocity for Aquarius halo Aq-A-5 run 
starting at redshift, z s , with Zel'dovich and 21pt resimulation ini- 
tial conditions. The results from the 21pt initial conditions are 
less sensitive to the start redshift. The cumulative function for 
the Zel'dovich initial conditions have been shifted upwards by 
one dex for clarity. The curve obtained by plotting the Zel'dovich 
z s = 127 with the shift lies within or close to the curves from the 
21pt initial conditions. 



In Figure [5] we show the circular velocity as a function 
of radius for Aq-A-5 halo at redshift zero for all eight simula- 
tions. The radius is taken from th e position of the pot ential 
minimum identified by SUBFIND jSpringel et al.ll200ll) . and 
the circular velocity is computed assuming spherical symme- 
try for the mass distribution about the potential minimum. 
Comparing the circular velocities for the Zel'dovich initial 
initial conditions amongst themselves, and likewise the 21pt 
circular velocities (displaced vertically by 30 kms^ 1 for clar- 
ity) it is first evident that the circular velocity curves are 
almost all rather alike, with the simulation starting from 
z a = 15 Zel'dovich initial conditions being the outlier. For 
these initial conditions the peak of the circular velocity at 
z = occurs almost 50% further from the centre than in 
the other simulations. Comparing the curves it is clear that 
the circular velocity shows a much weaker dependence on 
the start redshift with the 21pt initial conditions than with 
Zel'dovich initial conditions. 

In Figure HI] we plot the cumulative subhalo abundance 
as a function of the maximum subhalo circular velocity for 
the eight simulations, shifting the Zel'dovich curves upward 
by one dex for clarity. Again, the subhalo circular velocity 
functions show a much greater consistency at z=0 amongst 
the simulations starting from 21pt initial conditions com- 
pared with Zel'dovich initial conditions. 

Finally in this Subsection, we look at the sensitivity of 
the position of subhaloes inside the Aq-A-5 halo, as deter- 
mined by SUBFIND, to the start redshift and kind of initial 
conditions. A comparison of the relative subhalo positions 
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Figure 7. Sensitivity of the subhalo positions to the start red- 
shift and type of initial condition. Each circle above the diagonal 
represents the result of a comparison, made at z=0, between two 
resimulations out of a total of eight resimulations of halo Aq-A-5. 
The type of initial condition and the start redshift, z s , are given 
in the figure legend. The radius of each circle is equal to the ge- 
ometric mean distance between a sample of matched subhalos, 
common to all eight resimulations. For scale the circle in the leg- 
end has a radius of 192kpc//i which corresponds to the geometric 
mean separation, averaged over all eight resimulations, between 
randomly chosen pairs of subhalos within the same subhalo. The 
21pt initial conditions are significantly less sensitive to the start 
redshift than Zel'dovich initial conditions. 

for the Aq- A halo resimulated at different resolutions was 
presented in lSpringel et a l. (2008). For these simulations it 
proved possible to identify corresponding substructures in 
the different simulations and compare their relative posi- 
tions. All of these simulations were started from the same 
start redshift with Zel'dovich initial conditions so the details 
of the comparison are not directly relevant, nonetheless the 
fact that differences of similar magnitude to those found here 
occur, implies that there are other numerical issues which 
also affect the positions of substructures significantly. 

We can likewise match substructures between our sim- 
ulations, and this is made easier by the fact that there is a 
one to one correspondence between the dark matter parti- 
cles in the different resimulations based on their Lagrangian 
coordinates. To make the a fair comparison between all 28 
pairs of simulations we first found a subset of matched sub- 
halos common to all eight resimulations. This sample has 
only 50 members. This compares to between 102 and 180 
matches between the different pairs of simulations. The fail- 
ure to match can arise from two boundary effects: a subhalo 
in one simulation may be a free halo in the other; a 'subhalo' 
may be below the limit of 32 particles in one of the simula- 
tions. Almost half of the subhalos found have less than 64 
particles. 

The subhalo positions in two different simulations are 
compared, after first subtracting the position of the main 



subhalo. There are many ways to compare the differences 
found. We have selected to show the geometric mean of the 
distance between matched subhalos between pairs of simu- 
lations in Figure [7) Comparisons based on other quantities 
such as the mean or median separation give qualitatively 
similar results. 

Again we find, from Figure [7] that the simulations 
started from 21pt initial conditions show much less varia- 
tion in position on the start redshift than Zel'dovich initial 
conditions. For the simulation with Zel'dovich initial con- 
ditions starting from z s — 15, the location of the subhalos 
is very poorly correlated with any other of the simulations, 
although still better than the geometric mean separation of 
all non-identical pairs of subhalos in any given halo. Start- 
ing the Zel'dovich initial conditions at z s — 127, the starting 
redshift for the Aquarius halos, does give subhalo positions 
which match the 21pt resimulation subhalo positions fairly 
well. We checked that the analogous plot produced using 
the larger matched subhalos between each pair, rather than 
the subsample of subhalos common to all eight simulations, 
gives a very similar plot to the plot shown. 

In conclusion, for a resimulation of a Milky- Way like 
dark matter halo it is clear that the 2 = properties of the 
halo do depend somewhat on the start redshift, but that 
dependence is much reduced with 21pt resimulation initial 
conditions. The evolution of the halo to z = is highly non- 
linear, but the results are qualitatively in agreement with the 
known behaviour of 21pt and Zel'dovich init ial conditions in 
the quasi-linear regime (|Scoccimarrol Il998h where this can 
be understood from the behaviour of transients which are 
both smaller and decay faster (by a factor 1 + z) for 21pt 
initial conditions than for Zel'dovich initial conditions. In 
quantitative terms the differences amongst the circular ve- 
locity curves, the cumulative subhalo abundances and the 
relative positions of substructures between the different sim- 
ulations started from 21pt initial conditions, is larger than 
naively expect from the behaviour of transients. This is most 
likely because there are other numerical factors which lead 
to much of the differences seen in the various simulations 
starting from 21pt initial conditions. 

To show this is indeed the case we repeated one res- 
imulation using initial conditions that were mirror reflected 
about the centre of the periodic volume. Running these oth- 
erwise identical simulations leads, after reflecting the final 
conditions, to differences in the positions of the subhalos. 
The numerical process of mirroring is not a perfect sym- 
metry operation when operated on double precision float- 
ing point numbers so in effect the comparison shows that 
the effects of rounding error in the positions can propagate 
to much larger scales. After matching 189 subhalos for the 
z=0 halos we find the geometric separation of the subhalos 
between these reflected initial conditions is 12 kpc/h. This 
can be compared to 26 kpc//i difference between the 175 
matched subhalos between the 21pt resimulations starting 
at redshifts 127 and 63. This kind of sensitivity to the ini- 
tial conditions is expected given that halo formation is a 
strongly non-linear processes. It is interesting however, even 
though these simulations of halo formation show this kind 
of sensitivity to small changes in the initial conditions, that 
the effects of starting from 21pt or Zel'dovich initial condi- 
tions are still apparent at redshift zero for many properties 
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including for the circular velocity curve, the subhalo velocity 
functions and subhalo positions. 



5 SUMMARY AND DISCUSSION 

We have implemented a new method to make second or- 
der Lagrangian perturbation theory (21pt) initial conditions 
which is well suited to creating resimulation initial condi- 
tions. We have tested the method for an analytic spherically 
symmetric perturbation with features ranging over a factor 
of 5000 in linear scale or eleven orders of magnitude in mass, 
and can reproduce the analytic solution to a fractional ac- 
curacy of better than one percent over this range of scales 
in radius measured from the symmetry centre. 

Applying the new method we have recreated the ini- 
tial conditions for th e formation of a high redshift halo 
from Gao et all ||2005f ) and the a Milky- way mass halo from 
ISpringel et alj (|2008l ). We find from studying the properties 
of the final halos, that the final conditions show much less 
sensitivity to the start redshift when using 21pt initial con- 
ditions than when started from Zel'dovich initial conditions. 

For didactic purposes, we have calculated the effect of 
using Zel'dovich and 21pt initial conditions for a spherical 
top-hat collapse. In this simple case, the epoch of collapse 
is delayed by a timing error which depends primarily on the 
number of expansions between the epoch of the initial condi- 
tions and the collapse time. This timing error, for fixed start- 
ing redshift, grows with collapse time for Zel'dovich initial 
conditions, but decreases with collapse time for 21pt initial 
conditions. Applying the top-hat timing errors as a correc- 
tion to the radius of shells containing fixed amounts of mass 
as a function of time has the effect of bring ing the Zel ' dovich 
and 21pt resimulations of the high redshift iGao et all (|2005T ) 
halo into near coincidence. 

We find that for the Aq-A-5 halo the choi ce of start red- 
shift f or the Zel'dovich initial conditions in ISpringel et al.l 
(2008) was sufficiently high at z s = 127 to make little differ- 
ence to the properties of the halo at z = when compared to 
runs starting from 21pt initial conditions. For lower starting 
redshifts (z s — 63, 31, 15) the positions of substructures at 
z = do become sensit ive to the precise start redshift. 

Kncbe et alj (|2009l ) have also looked at the properties 
of halos, for the mass range (10 10 - 10 13 h~ 1 M Q ) at redshift 
zero, for simulations starting from Zel'dovich and 21pt ini- 
tial co nditions (generated using the code by ICrocce et al.1 
(2006)). In their paper they simulated a cosmological vol- 
ume which yielded many small halos and about ten halos 
with 150000 or more particles within the virial radius. In 
their study they looked at the distribution of halo proper- 
ties including the spin, triaxiality and concentration. They 
also matched halos in a similar way to the way we have 
matched subhalos and looked at the ratios of masses, tri- 
axialities, spins and concentrations of the matching halos. 
The authors concluded that any actual trends with start 
redshift or type of initial condition by redshift zero are cer- 
tainly small (for start redshifts of 25 - 150). It is not possible 
to directly compare our results with theirs, but the trends 
in the halo properties we have observed in halo Aq-A-5 at 
redshift zero are indeed small, and not obviously in disagree- 
ment with their results for populations of halos. 

For many purposes, the differences which arise from us- 



ing Zel'dovich or 21pt initial conditions may be small when 
compared to the uncertainties which arise in the modelling 
of more complex physical processes. From this point of view 
the main advantage of using 21pt initial conditions is that 
they allow the simulations to start later thus saving com- 
puter time. The issue of exactly what redshift one can start 
depends on the scientific problem, but the spherical collapse 
model, discussed in Subsection 12.21 gives a quantitative way 
to compare the start redshift for Zel'dovich or 21pt initial 
conditions. However it remains the case that it is advis- 
able to test the sensitivity of simulation results to the start 
redshi ft by direct simulation as for example in lPower et al.l 
(2003) for resimulations of galactic dark matter halos. 



ACKNOWLEDGEMENTS 

I would particularly like to thank Volker Springel for im- 
plementing part of the method described in this paper in 
his P-Gadget3 code and for reading a draft of this paper. I 
would also like to thank Roman Scoccimarro for clarifying 
an issue connected with the second-order overdensity field 
early o n in this project, Gao Liang, for helping with remak- 
ing the lGao et al.l l|2005h initial conditions and Carlos Frenk 
for comments on the first draft. The simulations in this pa- 
per were run on COSMA and I am grateful to Lydia Heck 
for computer support. This work was supported by STFC 
rolling grant ST/F002289/1. 



REFERENCES 

Aarseth S. J., Turner E. L., Gott III J. R., 1979, ApJ, 228, 
664 

Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, 

ApJ, 304, 15 
Bertschinger E., 2001, ApJS, 137, 1 

Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, 

A&A, 296, 575 
Buchert T., 1994, MNRAS, 267, 811 

Buchert T., Melott A. L., Weiss A. G., 1994, A&A, 288, 
349 

Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 
369 

Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, 
ApJ, 292, 371 

Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, 
ApJS, 57, 241 

Evrard A. E., MacFarland T. J., Couchman H. M. P., Col- 
berg J. M., Yoshida N., White S. D. M., Jenkins A., Frenk 
C. S., Pearce F. R., Peacock J. A., Thomas P. A., 2002, 
ApJ, 573, 7 

Gao L., White S. D. M., Jenkins A., Frenk C. S., Springel 
V., 2005, MNRAS, 363, 379 

Ghigna S., Moore B., Governato F., Lake C, Quinn T., 
Stadel J., 1998, MNRAS, 300, 146 

Hockney R. W., Eastwood J. W., 1988, Computer simula- 
tion using particles. Bristol: Hilger, 1988 

Hoffman Y., Ribak E., 1991, ApJ, 380, L5 

Joyce M., Marcos B., Gabrielli A., Baertschiger T., Sylos 
Labini F., 2005, Physical Review Letters, 95, 011304 



© 0000 RAS, MNRAS 000, 000-000 



14 A. Jenkins 



Katz N., Quinn T., Bertschinger E., Gelb J. M., 1994, MN- 
RAS, 270, L71 

Knebe A., Wagner C, Knollmann S., Diekershoff T., 

Krause F., 2009, ApJ, 698, 266 
Moore B., Ghigna S., Governato F., Lake G., Quinn T., 

Stadel J., Tozzi P., 1999, ApJ, 524, L19 
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 

563 

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 
493 

Navarro J. F., Hayashi E., Power C, Jenkins A. R., Frenk 
C. S., White S. D. M., Springel V., Stadel J, Quinn T. R., 
2004, MNRAS, 349, 1039 
Navarro J. F., White S. D. M., 1994, MNRAS, 267, 401 
Peebles P. J. E., 1980, The large-scale structure of the uni- 
verse 

Peebles P. J. E., 1982, ApJ, 263, LI 
Pen U.-L., 1997, ApJ, 490, L127 

Power C, Navarro J. F., Jenkins A., Frenk C. S., White 
S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 
338, 14 

Press W. H., Schechter P., 1974, ApJ, 187, 425 

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery 

B. P., 1992, Numerical recipes in FORTRAN. The art of 

scientific computing 
Prunet S., Pichon C, Aubert D., Pogosyan D., Teyssier R., 

Gottloeber S., 2008, ApJS, 178, 179 
Reed D. S., Bower R., Frenk C. S., Gao L., Jenkins A., 

Theuns T., White S. D. M., 2005, MNRAS, 363, 393 
Salmon J., 1996, ApJ, 460, 59 
Scoccimarro R., 1998, MNRAS, 299, 1097 
Sirko E., 2005, ApJ, 634, 728 
Springel V., 2005, MNRAS, 364, 1105 

Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins 
A., Helmi A., Navarro J. F., Frenk C. S., White S. D. M., 
2008, MNRAS, 391, 1685 

Springel V., White S. D. M., Jenkins A., Frenk C. S., 
Yoshida N., Gao L., Navarro J., Thacker R., Croton D., 
Helly J., Peacock J. A., Cole S., Thomas P., Couchman 
H., Evrard A., Colberg J., Pearce F., 2005, Nature, 435, 
629 

Springel V., White S. D. M., Tormen G., Kauffmann G., 

2001, MNRAS, 328, 726 
Stadel J., Potter D., Moore B., Diemand J., Madau P., 

Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, pp L280- 
White S. D. M., 1996, in Schaeffer R., Silk J., Spiro M., 

Zinn- Justin J., eds, Cosmology and Large Scale Structure 

Formation and Evolution of Galaxies, pp 349- 
Zel'dovich Y. B., 1970, A&A, 5, 84 



6 APPENDIX: METHODS FOR 

INTERPOLATION OF DISPLACEMENT 
FIELDS 

As mentioned in the main text, it is necessary to use interpo- 
lation of the displacement fields when making resimulation 
initial conditions in order to find the displacements at the 
particle positions. This need arises for example because the 
interparticle spacing in the high resolution region of the res- 
imulation initial conditions is typically many times smaller 



than the Fourier grid used to generate the longest wave- 
length displacements over the entire periodic simulation vol- 
ume. Another situation where interpolation is needed, even 
for cosmological initial conditions, occur s when the p article 
load is based on a glass-like distribution (|Whitell 19961 ) which 
means the particles cannot coincide with the grid points. 

In this Appendix we will examine two interpolation 
methods which are in use and compare them. We will sim- 
plify the discussion to the case in the limit of the inter- 
particle separation going to zero creating a continuous and 
uniform density medium. We assume that the displacement 
field is built of a finite number of plane waves and that 
the displacements are known precisely on a regular cubic 
grid of points. The purpose of the interpolation is to recon- 
struct the displacement field everywhere between these grid 
points. The success of a particular interpolation scheme can 
be judged by comparing the density field created by apply- 
ing the displacements given by interpolated displacement 
field with the density field generated from the true displace- 
ment field. We will assume the displacements are small so 
there is no shell crossing and further consider only the linear 
dependence of the density field on the displacements. The 
interpolation schemes we consider use linear combinations 
of the values known at the grid points, so it is sufficient to 
characterise the interpolation scheme by applying it to plane 
waves. 

We will compare the performance of two interpolation 
sche mes: the first, trilinear interpo lation or Cloud-in-Cell or 
CIC. iHocknev fc Eastwoodl |l988), has been used for mak- 
ing initial conditions over many years for ex ample in the 
Hubble Volume simu lations (I Evrard et al|2002l ) and the Mil- 
lennium simulation (|Springel et al.l 20051) as w ell as resim- 
ulation initial conditions e. g. lGao et alffeoOSl ). The second 
method, not previously documented, which we call QI, for 
quadratic interpolation, has been used more rece ntly for the 
Aquarius initial conditions (|Springel et al.ll2008l ). 

In principle, because the displacement fields for resim- 
ulations are generated by Fourier methods are bandwidth 
limited, the interpolation of the displacement field to the 
particle positions can be done pe rfectly using the W hittaker- 
Shannon interpolation formula (|Press et al.lll992l ). This is 
impractical because the sum must be taken over all grid 
points. The only exception to this is the trivial case where 
the particles lie on the Fourier grid points themselves, in 
which case no interpolation is needed. For cosmological ini- 
tial conditions, unlike resimulation initial conditions, this 
can easily be arranged, as for e xample in the 21pt initial 
condition code made available bv lCrocce et alj ([2006). 

Before discussing particular interpolation formulae we 
first define some convenient notation. Suppose the compu- 
tational domain is a periodic cube of side length L and that 
we define a cubic grid with p 3 cubic cells of side length A 
such that L — pA. A general point in the volume has co- 
ordinates (x, y, z) where all three spatial coordinates have a 
range < x,y, z < L. We will consider the interpolation of a 
scalar field A(x, y, z) and we will assume A a bandwidth lim- 
ited function, in such a way that its value can be predicted 
exactly anywhere within the domain using the Whittaker- 
Shannon interpolation formula and summing over the p 3 grid 
points. 

We first consider trilinear interpolation. The interpo- 
lated field, Atli, is given by the expression: 
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A TLI = W(^-h)W(^-h)W(^-h)A, (24) 

where the (ii, h, h) represents an integer triplet and the sum 
is over all p 3 grid points. For TLI the kernel W is defined 
as: 



W{a) = 



1 - \a\ if \a\ < 1 







if \a\ > 1. 



(25) 



TLI interpolation is exact for a function which is a sum 
of terms x ai y a2 z°" 3 where < on < 1. The compact form 
of W means that to evaluate the interpolated field at the 
position of a particle requires summing only eight terms 
corresponding to the nearest grid points. The resulting field, 
Atu, is continuous everywhere, which is a necessary require- 
ment when interpolating a component of the displacement 
field, otherwise regions of zero or double density can be cre- 
ated along the interfaces marking the discontinuities. When 
TLI is applied to the components of the displacement field, 
the resulting density field has discontinuities. 

We will now consider the effect on a plane wave. The 
'interpolated' wave will no longer be a simple plane wave 
but a superposition of the 'fundamental' plane wave (with a 
reduced amplitude) and the harmonics of the fundamental 
wave. The same is true of the density field generated by 
applying the interpolated displacements (assumed normal 
to the plane wave and of small amplitude) to a uniform 
medium. The relative power in the fundamental wave and 
the harmonics will differ between the displacement field and 
the density field, with relatively more power in the density 
field harmonics. The ideal interpolation scheme minimises 
the power in these harmonics and maximises the power in 
the fundamental plane wave. 

Consider a plane wave mode in a periodic cube of side- 
length L : 

2ir 

A(ki, k 2 , k 3 ) = Bexp[i — (kix + kiy + k^z)], (26) 

Li 

where fei, k 2 , kz are integers and B is the amplitude of the 
wave. It is straightforward to show that applying TLI inter- 
polation to the plane wave gives a reduced amplitude, -Btli, 
for the fundamental plain wave given by the following ratio: 



Btia/B = 



sin u± 

Ul 



smit2 

U2 



u 3 



(27) 



where Ui = nki/p provided < Ui < n/2. 

This factor applies equally to the displacement field and 
the linear density field created by applying the interpolated 
displacement field to a uniform medium. For small values of 
Ui the ratio is close unity and deviates below unity quadrat- 
ically to lowest order. It is desirable to use as large a Fourier 
grid as memory allows for making resimulation initial con- 
ditions to minimise the interpolation error. 

The interpolation can be made more accurate if the 
derivatives of the field, A, are also known at the grid points. 
The derivatives of the displacement field need to be com- 
puted for 21pt initial conditions in any case and it was the 
fact that the derivatives need to be calculated anyway that 
motivated trying to find an improved interpolation scheme 
making use of them. 

It is straightforward, making use of the derivatives of 



the field to devise a formula which is exact for linear and 
quadratic functions of the coordinates. We will call it the 
quadratic scheme. The field ^4qi is given by: 

^4qi =j4tli+ 



(28) 
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where: 



o(l-o)/2 ifO<a<l 
V(a) = { -a(l + a)/2 if -1 < a < 
if \a\ > 1. 



(29) 



The QI interpolation is exact for any function which is a sum 
of terms x ai y a2 z a3 where < on < 2. The three additional 
terms in the QI formula vanish at the grid points so like 
TLI interpolation formula, the QI scheme is exact at the 
grid points. The corresponding effect of QI interpolation on 
a plane wave can be shown to be: 



Bqi/B = 




(30) 



where m = irki/p provided < m < 7r/2. Again this ratio is 
close to unity, when m is small, but the leading order devia- 
tion below unity for small Uj is now quartic - a considerable 
improvement. 

Figure [8] compares the two interpolation schemes act- 
ing on a plane wave along one of the principal coordinate 
axes. The upper curves show the ratio of the power in the 
density field created by applying the interpolated displace- 
ments with the true power for both TLI and QI interpolation 
schemes. The curves are obtained by setting U2 — 0; u$ — 
in equations (|27[) and (|30[) respectively and the x-axis is 
plotted up to the Nyquist wavenumber of the Fourier grid 
points where the displacement field is known exactly. Both 
interpolation schemes introduce power into the harmonics 
of the plane wave. For a 1-d plane wave the TLI scheme this 
harmonic power can be shown to be: 



(31) 



while for the QI scheme the harmonic power is given by: 



-Pha 



Ptr 



u 



2 sin u sin 2u 



■2a 



(32) 

where u — kn/p provided < u < tv/2. The lower two curves 
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Figure 8. The curves show the effects of the TLI and QI interpo- 
lation schemes when applied to a 1-d on a plane wave, represent- 
ing a displacement field normal to the wave, of wavenumbcr k, 
given a sampling rate of the plane wave of wavenumber ^Nyquist • 
The upper curves show the power remaining in resulting linear 
density created by displacing a uniform medium, with the plane 
wave displacements given by interpolation, while the lower curves 
show the power introduced into the the harmonics of the plane 
wave by the process of interpolation. 



show the amount of power generated in the harmonics for 
the TLI and QI interpolation schemes. The QI scheme is 
evidently superior to TLI on both measures. 

Both the TLI and QI interpolation schemes have cubic 
symmetry, and therefore generate some anisotropies. How- 
ever in both schemes these anisotropies are small and enter 
first at quartic order in the magnitude of the wavevector. 
The fractional RMS derivation in the power in the funda- 
mental wave about the average over all directions less than 
1 percent even at the Nyquist wavelength for both TLI and 
QI, with QI being the best by a factor three. 

For a fixed grid size implementing the QI formula is a 
factor of four more expensive than TLI as the three deriva- 
tives of A need to be computed in addition to A itself. For 
initial conditions, this factor is reduced to three, because 
the quantities of interest are gradients of a potential which 
reduces the number of independent terms that have to be 
evaluated. The QI scheme however is more attractive after 
taking account of its improved accuracy and the fact that 
the computational effort scales approximately as the cube 
of the dimension of the Fourier grid. For limited memory or 
limited CPU the QI scheme is to be preferred to TLI. 

The QI interpolation scheme was used for making all 
the Zel'dovich and 21pt initial conditions used in this paper. 
For the derivatives of the Zel'dovich displacements, which 
are required to compute the second-order overdensity, using 
the QI interpolation scheme entails calculating the second 
order derivatives of the Zel'dovich displacements. 
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